/************************** NCHYPPRO.CPP ********************** 2002-10-20 AF *
*
* Calculation of univariate and multivariate noncentral hypergeometric
* probability distribution.
*
* This file contains source code for the class NonCentralHypergeometricC
* and MultiNonCentralHypergeometricMoments defined in stocc.h.
*
* Documentation:
* ==============
* The file stocc.h contains class definitions.
* The file stocc.htm contains further instructions.
* The file nchyp.pdf, available from www.agner.org/random/theory
* describes the theory of the calculation methods.
*
*******************************************************************************/

#include <math.h>      // math functions
#include <string.h>    // memcpy function
#include "stocc.hpp"   // class definition
#include "erfres.cpp"  // table of error function residues for Laplace's method

/***********************************************************************
         Shared functions
***********************************************************************/

double pow2_1( double q, double * y0 = 0 ) {
   // calculate 2^q and (1-2^q) without loss of precision.
   // return value is (1-2^q). 2^q is returned in *y0
   double y, y1, y2, qn, i, ifac;
   q *= M_LN2;
   if ( fabs( q ) > 0.1 ) {
      y = exp( q );
      y1 = 1. - y;
   }
   else { // expand 1-e^q = -summa(q^n/n!) to avoid loss of precision
      y1 = 0;
      qn = i = ifac = 1;
      do {
         y2 = y1;
         qn *= q;
         ifac *= i++;
         y1 -= qn / ifac;
      } while ( y1 != y2 );
      y = 1.-y1;
   }
   if ( y0 ) *y0 = y;
   return y1;
}

double pow1pow( double q, double x ) {
   // calculate ln((1-e^q)^x) without loss of precision
   // Uses various Taylor expansions to avoid loss of precision
   double y, y1, y2, z1, z2, qn, i, ifac;

   if ( fabs( q ) > 0.1 ) {
      y = exp( q );
      y1 = 1. - y;
   }
   else { // expand 1-e^q = -summa(q^n/n!) to avoid loss of precision
      y1 = 0;
      qn = i = ifac = 1;
      do {
         y2 = y1;
         qn *= q;
         ifac *= i++;
         y1 -= qn / ifac;
      } while ( y1 != y2 );
      y = 1. - y1;
   }
   if ( y > 0.1 ) { // (1-y)^x calculated without problem
      return x * log( y1 );
   }
   else { // expand ln(1-y) = -summa(y^n/n)
      y1 = i = 1.;
      z1 = 0;
      do {
         z2 = z1;
         y1 *= y;
         z1 -= y1 / i++;
      } while ( z1 != z2 );
      return x * z1;
   }
}

double log1mx( double x, double x1 ) {
   // calculate log(1-x) without loss of precision when x is small
   // parameter x1 must be = 1-x
   if ( x > 0.03 ) {
      return log( x1 );
   }
   else { // expand ln(1-x) = -summa(x^n/n)
      double y, z1, z2, i;
      y = i = 1.;
      z1 = 0;
      do {
         z2 = z1;
         y *= x;
         z1 -= y / i++;
      } while ( z1 != z2 );
      return z1;
   }
}

double LnFacr( double x ) {
   // log factorial of non-integer x
   long int ix = ( long int )( x );
   if ( x == ix ) return LnFac( ix ); // x is integer
   double r, r2, D = 1., f;
   static const double
   C0 =  0.918938533204672722,   // ln(sqrt(2*pi))
   C1 =  1./12.,
   C3 = -1./360.,
   C5 =  1./1260.,
   C7 = -1./1680.;
   if ( x < 6. ) {
      if ( x == 0 || x == 1 ) return 0;
      while ( x < 6 ) D *= ++x;
   }
   r  = 1. / x;
   r2 = r*r;
   f = ( x + 0.5 )*log( x ) - x + C0 + r*( C1 + r2*( C3 + r2*( C5 + r2*C7 ) ) );
   if ( D != 1. ) f -= log( D );
   return f;
}

double FallingFactorial( double a, double b ) {
   // calculates ln(a*(a-1)*(a-2)* ... * (a-b+1))

   if ( b < 30 && int( b ) == b && a < 1E10 ) {
      // direct calculation
      double f = 1.;
      for ( int i=0; i<b; i++ ) f *= a--;
      return log( f );
   }

   if ( a > 100.*b && b > 1. ) {
      // combine Stirling formulas for a and (a-b) to avoid loss of precision
      double ar = 1./a;
      double cr = 1./( a-b );
      // calculate -log(1-b/a) by Taylor expansion
      double s = 0., lasts, n = 1., ba = b*ar, f = ba;
      do {
         lasts = s;
         s += f/n;
         f *= ba;
         n++;
      } while ( s != lasts );
      return ( a+0.5 )*s + b*log( a-b ) - b + ( 1./12. )*( ar-cr )
             //- (1./360.)*(ar*ar*ar-cr*cr*cr)
             ;
   }

   // use LnFacr function
   return LnFacr( a )-LnFacr( a-b );
}


/***********************************************************************
             Methods for class NonCentralHypergeometricC
***********************************************************************/

NonCentralHypergeometricC::NonCentralHypergeometricC( long int n_, long int m_, long int N_, double odds, double accuracy_ ) {
   // constructor
   n = n_;
   m = m_;
   N = N_;
   omega = odds;
   accuracy = accuracy_; // set parameters
   xmin = m + n - N;
   if ( xmin < 0 ) xmin = 0; // calculate xmin and xmax
   xmax = n;
   if ( xmax > m ) xmax = m;
   ParametersChanged = 1;
   r = 1.;
}     // initialize


double NonCentralHypergeometricC::mean( void ) {
   // find approximate mean
   int iter;                            // number of iterations
   double a, b;                         // temporaries in calculation of first guess
   double mean, mean1;                  // iteration value of mean
   double m1r, m2r;                     // 1/m, 1/m2
   double e1, e2;                       // temporaries
   double g;                            // function to find root of
   double gd;                           // derivative of g
   double omegar;                       // 1/omega

   if ( omega == 1. ) { // simple hypergeometric
      return double( m )*n/N;
   }

   // calculate Cornfield mean of extended hypergeometric distribution as first guess
   a = ( m+n )*omega + ( N-m-n );
   b = a*a - 4.*omega*( omega-1. )*m*n;
   b = b > 0. ? sqrt( b ) : 0.;
   mean = ( a-b )/( 2.*( omega-1. ) );
   if ( mean < xmin ) mean = xmin;
   if ( mean > xmax ) mean = xmax;

   iter = 0;
   m1r = 1./m;
   m2r = 1./( N-m );
   omegar = 1./omega;

   if ( omega > 1. ) {
      do { // Newton Raphson iteration
         mean1 = mean;
         e1 = 1.-( n-mean )*m2r;
         if ( e1 < 1E-14 ) {
            e2 = 0.;
         } // avoid underflow
         else {
            e2 = pow( e1,omega-1. );
         }
         g = e2*e1 + ( mean-m )*m1r;
         gd = e2*omega*m2r + m1r;
         mean -= g / gd;
         if ( mean < xmin ) mean = xmin;
         if ( mean > xmax ) mean = xmax;
         if ( ++iter > 40 ) {
            FatalError( "Search for mean failed" );
         }
      } while ( fabs( mean1 - mean ) > 2E-6 );
   }
   else { // omega < 1
      do { // Newton Raphson iteration
         mean1 = mean;
         e1 = 1.-mean*m1r;
         if ( e1 < 1E-14 ) {
            e2 = 0.;
         } // avoid underflow
         else {
            e2 = pow( e1,omegar-1. );
         }
         g = 1.-( n-mean )*m2r-e2*e1;
         gd = e2*omegar*m1r + m2r;
         mean -= g / gd;
         if ( mean < xmin ) mean = xmin;
         if ( mean > xmax ) mean = xmax;
         if ( ++iter > 40 ) {
            FatalError( "Search for mean failed" );
         }
      } while ( fabs( mean1 - mean ) > 2E-6 );
   }
   return mean;
}


double NonCentralHypergeometricC::moments( double * mean_, double * var_ ) {
   // calculate exact mean and variance
   // return value = sum of f(x), expected = 1.
   double y, sy=0, sxy=0, sxxy=0, me1;
   long int x, xm, x1;
   const float accuracy = 1E-10f;  // accuracy of calculation
   xm = ( long int )mean(); // approximation to mean
   for ( x=xm; x<=xmax; x++ ) {
      y = probability( x );
      x1 = x - xm;  // subtract approximate mean to avoid loss of precision in sums
      sy += y;
      sxy += x1 * y;
      sxxy += x1 * x1 * y;
      if ( y < accuracy ) break;
   }
   for ( x=xm-1; x>=xmin; x-- ) {
      y = probability( x );
      x1 = x - xm;  // subtract approximate mean to avoid loss of precision in sums
      sy += y;
      sxy += x1 * y;
      sxxy += x1 * x1 * y;
      if ( y < accuracy ) break;
   }

   me1 = sxy / sy;
   *mean_ = me1 + xm;
   y = sxxy / sy - me1 * me1;
   if ( y < 0 ) y=0;
   *var_ = y;
   return sy;
}


double NonCentralHypergeometricC::lnbico() {
   // natural log of binomial coefficients.
   // returns lambda = log(m!*x!/(m-x)!*m2!*x2!/(m2-x2)!)
   long int x2 = n-x, m2 = N-m;
   if ( ParametersChanged ) {
      mFac = LnFac( m ) + LnFac( m2 );
      xLast = -99;
      ParametersChanged = 0;
   }
   if ( m < FAK_LEN && m2 < FAK_LEN )  goto DEFLT;
   switch ( x - xLast ) {
   case 0: // x unchanged
      break;
   case 1: // x incremented. calculate from previous value
      xFac += log ( double( x ) * ( m2-x2 ) / ( double( x2+1 )*( m-x+1 ) ) );
      break;
   case -1: // x decremented. calculate from previous value
      xFac += log ( double( x2 ) * ( m-x ) / ( double( x+1 )*( m2-x2+1 ) ) );
      break;
   default:
DEFLT: // calculate all
      xFac = LnFac( x ) + LnFac( x2 ) + LnFac( m-x ) + LnFac( m2-x2 );
   }
   xLast = x;
   return bico = mFac - xFac;
}

void NonCentralHypergeometricC::findpars() {
   // calculate d, E, r, w

   // find r to center peak of integrand at 0.5
   double dd, d1, z, zd, rr, lastr, rrc, rt, r2, r21, a, b;
   double oo[2];
   double xx[2] = {x, n-x};
   int i, j = 0;
   if ( omega > 1. ) { // make both omegas <= 1 to avoid overflow
      oo[0] = 1.;
      oo[1] = 1./omega;
   } else {
      oo[0] = omega;
      oo[1] = 1.;
   }

   dd = oo[0]*( m-x ) + oo[1]*( N-m-xx[1] );
   d1 = 1./dd;
   E = ( oo[0]*m + oo[1]*( N-m ) ) * d1;
   rr = r;
   if ( rr <= d1 ) rr = 1.2*d1;         // initial guess
   // Newton-Raphson iteration to find r
   do {
      lastr = rr;
      rrc = 1. / rr;
      z = dd - rrc;
      zd = rrc * rrc;
      for ( i=0; i<2; i++ ) {
         rt = rr * oo[i];
         if ( rt < 100. ) {               // avoid overflow if rt big
            r21 = pow2_1( rt, &r2 );       // r2=2^r, r21=1.-2^r
            a = oo[i] / r21;               // omegai/(1.-2^r)
            b = xx[i] * a;                 // x*omegai/(1.-2^r)
            z  += b;
            zd += b * a * M_LN2 * r2;
         }
      }
      if ( zd == 0 ) FatalError( "can't find r" );
      rr -= z / zd;
      if ( rr <= d1 ) rr = lastr * 0.125 + d1*0.875;
      if ( ++j == 70 ) FatalError( "convergence problem searching for r" );
   } while ( fabs( rr-lastr ) > rr * 1.E-6 );
   if ( omega > 1 ) {
      dd *= omega;
      rr *= oo[1];
   }
   r = rr;
   rd = rr * dd;

   // find peak width
   double ro, k1, k2;
   ro = r * omega;
   if ( ro < 300 ) {                    // avoid overflow
      k1 = pow2_1( ro );
      k1 = -1. / k1;
      k1 = omega*omega*( k1+k1*k1 );
   } else k1 = 0.;
   if ( r < 300 ) {                     // avoid overflow
      k2 = pow2_1( r );
      k2 = -1. / k2;
      k2 = ( k2+k2*k2 );
   } else k2 = 0.;
   phi2d = -4.*r*r*( x*k1 + ( n-x )*k2 );
   if ( phi2d > 0. ) FatalError( "peak width undefined in function findpars" );
   w = 1./sqrt( -phi2d );
}


/***********************************************************************
         calculation methods in class NonCentralHypergeometricC
***********************************************************************/

double NonCentralHypergeometricC::recursive() {
   // recursive calculation
   // noncentral hypergeometric distribution by recursion formula
   // Approximate by ignoring probabilities < accuracy and minimize storage requirement
   const int BUFSIZE = 512;            // buffer size
   double accuracya;                   // absolute accuracy
   double p[BUFSIZE+2];                // probabilities
   double * p1, * p2;                  // offset into p
   double mxo;                         // (m-x)*omega
   double Nmnx;                        // N-m-nu+x
   double y, y1;                       // save old p[x] before it is overwritten
   double d1, d2, dcom;                // divisors in probability formula
   long int xi, nu;                    // xi, nu = recursion values of x, n
   long int x1, x2;                    // xi_min, xi_max

   accuracya = 0.005 * accuracy;       // absolute accuracy
   p1 = p2 = p + 1;                    // make space for p1[-1]
   p1[-1] = 0.;
   p1[0]  = 1.;          // initialize for recursion
   x1 = x2 = 0;
   for ( nu=1; nu<=n; nu++ ) {
      if ( n - nu < x - x1 || p1[x1] < accuracya ) {
         x1++;               // increase lower limit when breakpoint passed or probability negligible
         p2--;
      }              // compensate buffer offset in order to reduce storage space
      if ( x2 < x && p1[x2] >= accuracya ) {
         x2++;
         y1 = 0.;
      }    // increase upper limit until x has been reached
      else {
         y1 = p1[x2];
      }
      if ( x1 > x2 ) return 0.;
      if ( p2+x2-p > BUFSIZE ) FatalError( "buffer overrun in function recursive" );

      mxo = ( m-x2 )*omega;
      Nmnx = N-m-nu+x2+1;
      for ( xi = x2; xi >= x1; xi-- ) {  // backwards loop
         d2 = mxo + Nmnx;
         mxo += omega;
         Nmnx--;
         d1 = mxo + Nmnx;
         dcom = 1. / ( d1 * d2 );         // save a division by making common divisor
         y  = p1[xi-1]*mxo*d2*dcom + y1*( Nmnx+1 )*d1*dcom;
         y1 = p1[xi-1];                   // (warning: pointer alias, can't swap instruction order)
         p2[xi] = y;
      }
      p1 = p2;
   }

   if ( x < x1 || x > x2 ) return 0.;
   return p1[x];
}


double NonCentralHypergeometricC::binoexpand() {
   // calculate by binomial expansion of integrand
   // only for x < 2 or n-x < 2 (not implemented for higher x because of loss of precision)
   long int x1, m1, m2;
   double o;
   if ( x > n/2 ) { // invert
      x1 = n-x;
      m1 = N-m;
      m2 = m;
      o = 1./omega;
   } else {
      x1 = x;
      m1 = m;
      m2 = N-m;
      o = omega;
   }

   if ( x1 == 0 ) {
      return exp( FallingFactorial( m2,n ) - FallingFactorial( m2+o*m1,n ) );
   }

   if ( x1 == 1 ) {
      double d, e, q, q0, q1;
      q = FallingFactorial( m2,n-1 );
      e = o*m1+m2;
      q1 = q - FallingFactorial( e,n );
      e -= o;
      q0 = q - FallingFactorial( e,n );
      d = e - ( n-1 );
      return m1*d*( exp( q0 ) - exp( q1 ) );
   }

   FatalError( "x > 1 not supported by function binoexpand" );
   return 0;
}


double NonCentralHypergeometricC::laplace() {
   // Laplace's method with narrow integration interval,
   // using error function residue table, defined in erfres.cpp
   // NOTE: This function can only be used when the integrand peak is narrow.
   // findpars() must be called before this function.

   const int COLORS = 2;         // number of colors
   const int MAXDEG = 40;        // arraysize
   int degree;                   // max expansion degree
   double accur;                 // stop expansion when terms below this threshold
   double omegai[COLORS] = {omega, 1.}; // weights for each color
   double xi[COLORS] = {x, n-x}; // number of each color sampled
   double f0;                    // factor outside integral
   double rho[COLORS];           // r*omegai
   double qi;                    // 2^(-rho)
   double qi1;                   // 1-qi
   double qq[COLORS];            // qi / qi1
   double eta[COLORS+1][MAXDEG+1]; // eta coefficients
   double phideri[MAXDEG+1];     // derivatives of phi
   double PSIderi[MAXDEG+1];     // derivatives of PSI
   double epsilon;               // factor for integration interval
   double epsilonmax;            // max possible epsilon
   double * erfres;              // pointer to table of error function expansion residues
   int erfreslen;                // length of table

   // variables in asymptotic summation
   static const double sqrtpi = 1.772453850905515993; // sqrt(pi)
   static const double sqrt8  = 2.828427124746190098; // sqrt(8)
   double qqpow;                 // qq^j
   double pow2k;                 // 2^k
   double mfac;                  // (-1)^(k-1)*(k-1)!
   double kfac;                  // k!
   double bino;                  // binomial coefficient
   double wr;                    // 1/w, w = integrand peak width
   double vr;                    // 1/v, v = integration interval
   double v2m2;                  // (2*v)^(-2)
   double v2mk1;                 // (2*v)^(-k-1)
   double gammal2;               // Gamma((k+1)/2);
   double s;                     // summation term
   double sum;                   // Taylor sum

   int i;                        // loop counter for color
   int j;                        // loop counter for derivative
   int k;                        // loop counter for expansion degree
   int ll;                       // k/2
   int converg=0;                // number of consequtive terms below accuracy


   // initialize
   for ( k=0; k<=2; k++ )  phideri[k] = PSIderi[k] = 0;

   // find rho[i], qq[i], first eta coefficients, and zero'th derivative of phi
   for ( i=0; i<COLORS; i++ ) {
      rho[i] = r * omegai[i];
      if ( rho[i] > 40. ) {
         qi=0.;
         qi1 = 1.;
      }  // avoid underflow
      else {
         qi1 = pow2_1( -rho[i], &qi );
      }   // qi=2^(-rho), qi1=1.-2^(-rho)
      qq[i] = qi / qi1;    // 2^(-r*omegai)/(1.-2^(-r*omegai))
      // peak = zero'th derivative
      phideri[0] += xi[i] * log1mx( qi, qi1 );
      // eta coefficients
      eta[i][0] = 0.;
      eta[i][1] = eta[i][2] = rho[i]*rho[i];
   }

   // r, rd, and w must be calculated by findpars()
   // zero'th derivative
   phideri[0] -= ( rd - 1. ) * M_LN2;
   // scaled factor outside integral
   f0 = rd * exp( phideri[0] + lnbico() );
   // calculate narrowed integration interval
   vr = sqrt8 * w;
   wr = 1. / w;
   phideri[2] = phi2d;
   epsilonmax = 0.5 * wr;
   if ( accuracy < 1.E-10 && epsilonmax >= 8. ) {
      epsilon = 8.;
      erfres = erfres80;
      erfreslen = sizeof( erfres80 )/sizeof( erfres80[0] );
   } else if ( accuracy < 1.E-8 && epsilonmax >= 7. ) {
      epsilon = 7.;
      erfres = erfres70;
      erfreslen = sizeof( erfres70 )/sizeof( erfres70[0] );
   } else if ( accuracy < 1.E-6 && epsilonmax >= 6. ) {
      epsilon = 6.;
      erfres = erfres60;
      erfreslen = sizeof( erfres60 )/sizeof( erfres60[0] );
   } else {
      epsilon = 5.;
      erfres = erfres50;
      erfreslen = sizeof( erfres50 )/sizeof( erfres50[0] );
   }
   if ( epsilon > epsilonmax ) FatalError( "Laplace failed. peak width too high" );
   degree = MAXDEG;
   if ( degree >= erfreslen*2 ) degree = erfreslen*2-2;

   // variables used in summation
   v2m2 = 0.25*vr*vr;        // (2*v)^(-2)

   // set up for starting loop at k=3
   PSIderi[0] = 1.;
   pow2k = 8.;
   mfac = 2.;
   kfac = 2.;
   sum = 0.5 * vr * sqrtpi * erfres[0];
   v2mk1 = 0.5*vr*v2m2*v2m2;
   gammal2 = sqrtpi * 0.75;
   accur = accuracy * sum;

   // summation loop
   for ( k=3; k <= degree; k++ ) {
      phideri[k] = 0.;

      // loop for all (2) colors
      for ( i=0; i<COLORS; i++ ) {
         eta[i][k] = 0.;
         // backward loop for all powers
         for ( j=k; j>0; j-- ) {
            // find coefficients recursively from previous coefficients
            eta[i][j]  =  eta[i][j]*( j*rho[i]-( k-2 ) ) +  eta[i][j-1]*rho[i]*( j-1 );
         }
         qqpow = 1.;
         // forward loop for all powers
         for ( j=1; j<=k; j++ ) {
            qqpow *= qq[i];   // qq^j
            // contribution to derivative
            phideri[k] += xi[i] * eta[i][j] * qqpow;
         }
      }

      // finish calculation of derivatives
      phideri[k] = -pow2k * phideri[k] + 2*( 1-k )*phideri[k-1];

      pow2k *= 2.;    // 2^k
      mfac *= -k;     // (-1)^(k-1)*(k-1)!

      // loop to calculate derivatives of PSI from derivatives of psi.
      // terms # 0, 1, 2, k-2, and k-1 are zero and not included in loop.
      // The j'th derivatives of psi are identical to the derivatives of phi for j>2, and
      // zero for j=1,2. Hence we are using phideri[j] for j>2 here.
      PSIderi[k] = phideri[k]; // this is term # k
      bino = 0.5*( k-1 )*( k-2 ); // binomial coefficient for term # 3
      for ( j=3; j < k-2; j++ ) { // loop for remaining nonzero terms (if k>5)
         PSIderi[k] += PSIderi[k-j] * phideri[j] * bino;
         bino *= double( k-j )/double( j );
      }

      kfac *= k;

      if ( ( k & 1 ) == 0 ) { // only for even k
         ll = k/2;
         s = ( PSIderi[k]/kfac )*v2mk1* gammal2 * erfres[ll];
         sum += s;

         // check for convergence of Taylor expansion
         if ( fabs( s ) < accur ) converg++;
         else converg = 0;
         if ( converg > 1 ) break;

         // update recursive expressions
         v2mk1 *= v2m2;
         gammal2 *= k/2 + 0.5;
      }
   }

   // multiply by terms outside integral
   return f0 * sum;
}


double NonCentralHypergeometricC::integrate() {
   // non-central hypergeometric distribution function
   // calculation by numerical integration with variable-length steps
   // NOTE: findpars() must be called before this function.
   double s;                            // result of integration step
   double sum;                          // integral
   double ta, tb;                       // subinterval for integration step

   lnbico();                            // compute log of binomial coefficients

   // choose method:
   if ( w < 0.02 || w < 0.1 && ( x==m || n-x==N-m ) && accuracy > 1E-6 ) {
      // normal method. Step length determined by peak width w
      double delta, s1;
      s1 = accuracy < 1E-9 ? 0.5 : 1.;
      delta = s1 * w;                    // integration steplength
      ta = 0.5 + 0.5 * delta;
      sum = integrate_step( 1.-ta, ta ); // first integration step around center peak
      do {
         tb = ta + delta;
         if ( tb > 1. ) tb = 1.;
         s  = integrate_step( ta, tb );   // integration step to the right of peak
         s += integrate_step( 1.-tb,1.-ta ); // integration step to the left of peak
         sum += s;
         if ( s < accuracy * sum ) break; // stop before interval finished if accuracy reached
         ta = tb;
         if ( tb > 0.5 + w ) delta *= 2.;
      }  // increase step length far from peak
      while ( tb < 1. );
   }

   else {
      // difficult situation. Step length determined by inflection points
      double t1, t2, tinf, delta, delta1;
      sum = 0.;
      // do left and right half of integration interval separately:
      for ( t1=0., t2=0.5; t1 < 1.; t1+=0.5, t2+=0.5 ) {
         // integrate from 0 to 0.5 or from 0.5 to 1
         tinf = search_inflect( t1, t2 ); // find inflection point
         delta = tinf - t1;
         if ( delta > t2 - tinf ) delta = t2 - tinf; // distance to nearest endpoint
         delta *= 1./7.;                  // 1/7 will give 3 steps to nearest endpoint
         if ( delta < 1E-4 ) delta = 1E-4;
         delta1 = delta;
         // integrate from tinf forwards to t2
         ta = tinf;
         do {
            tb = ta + delta1;
            if ( tb > t2 - 0.25*delta1 ) tb = t2; // last step of this subinterval
            s = integrate_step( ta, tb );  // integration step
            sum += s;
            delta1 *= 2;                   // double steplength
            if ( s < sum * 1E-4 ) delta1 *= 8.; // large step when s small
            ta = tb;
         } while ( tb < t2 );
         if ( tinf ) {
            // integrate from tinf backwards to t1
            tb = tinf;
            do {
               ta = tb - delta;
               if ( ta < t1 + 0.25*delta ) ta = t1; // last step of this subinterval
               s = integrate_step( ta, tb ); // integration step
               sum += s;
               delta *= 2;                  // double steplength
               if ( s < sum * 1E-4 ) delta *= 8.; // large step when s small
               tb = ta;
            } while ( ta > t1 );
         }
      }
   }

   return sum * rd;
}


double NonCentralHypergeometricC::integrate_step( double ta, double tb ) {
   // integration subprocedure used by integrate()
   // makes one integration step from ta to tb using Gauss-Legendre method.
   // result is scaled by multiplication with exp(bico)
   double ab, delta, tau, ltau, y, sum, taur, rdm1;
   int i;

   // define constants for Gauss-Legendre integration with IPOINTS points
#define IPOINTS  8  // number of points in each integration step

#if IPOINTS == 3
   static const double xval[3]    = {-.774596669241,0,0.774596668241};
   static const double weights[3] = {.5555555555555555,.88888888888888888,.55555555555555};
#elif IPOINTS == 4
   static const double xval[4]    = {-0.861136311594,-0.339981043585,0.339981043585,0.861136311594},
                                    static const double weights[4] = {0.347854845137,0.652145154863,0.652145154863,0.347854845137};
#elif IPOINTS == 5
   static const double xval[5]    = {-0.906179845939,-0.538469310106,0,0.538469310106,0.906179845939};
   static const double weights[5] = {0.236926885056,0.478628670499,0.568888888889,0.478628670499,0.236926885056};
#elif IPOINTS == 6
   static const double xval[6]    = {-0.932469514203,-0.661209386466,-0.238619186083,0.238619186083,0.661209386466,0.932469514203};
   static const double weights[6] = {0.171324492379,0.360761573048,0.467913934573,0.467913934573,0.360761573048,0.171324492379};
#elif IPOINTS == 8
   static const double xval[8]    = {-0.960289856498,-0.796666477414,-0.525532409916,-0.183434642496,0.183434642496,0.525532409916,0.796666477414,0.960289856498};
   static const double weights[8] = {0.10122853629,0.222381034453,0.313706645878,0.362683783378,0.362683783378,0.313706645878,0.222381034453,0.10122853629};
#elif IPOINTS == 12
   static const double xval[12]   = {-0.981560634247,-0.90411725637,-0.769902674194,-0.587317954287,-0.367831498998,-0.125233408511,0.125233408511,0.367831498998,0.587317954287,0.769902674194,0.90411725637,0.981560634247};
   static const double weights[12]= {0.0471753363866,0.106939325995,0.160078328543,0.203167426723,0.233492536538,0.249147045813,0.249147045813,0.233492536538,0.203167426723,0.160078328543,0.106939325995,0.0471753363866};
#elif IPOINTS == 16
   static const double xval[16]   = {-0.989400934992,-0.944575023073,-0.865631202388,-0.755404408355,-0.617876244403,-0.458016777657,-0.281603550779,-0.0950125098376,0.0950125098376,0.281603550779,0.458016777657,0.617876244403,0.755404408355,0.865631202388,0.944575023073,0.989400934992};
   static const double weights[16]= {0.027152459411,0.0622535239372,0.0951585116838,0.124628971256,0.149595988817,0.169156519395,0.182603415045,0.189450610455,0.189450610455,0.182603415045,0.169156519395,0.149595988817,0.124628971256,0.0951585116838,0.0622535239372,0.027152459411};
#else
#error // IPOINTS must be a value for which the tables are defined
#endif

   delta = 0.5 * ( tb - ta );
   ab = 0.5 * ( ta + tb );
   rdm1 = rd - 1.;
   sum = 0;

   for ( i = 0; i < IPOINTS; i++ ) {
      tau = ab + delta * xval[i];
      ltau = log( tau );
      taur = r * ltau;
      // possible loss of precision due to subtraction here:
      y = pow1pow( taur*omega,x ) + pow1pow( taur,n-x ) + rdm1*ltau + bico;
      if ( y > -50. ) sum += weights[i] * exp( y );
   }
   return delta * sum;
}


double NonCentralHypergeometricC::search_inflect( double t_from, double t_to ) {
   // search for an inflection point of the integrand PHI(t) in the interval
   // t_from < t < t_to
   const int COLORS = 2;                // number of colors
   double t, t1;                        // independent variable
   double rho[COLORS];                  // r*omega[i]
   double q;                            // t^rho[i] / (1-t^rho[i])
   double q1;                           // 1-t^rho[i]
   double xx[COLORS];                   // x[i]
   double zeta[COLORS][4][4];           // zeta[i,j,k] coefficients
   double phi[4];                       // derivatives of phi(t) = log PHI(t)
   double Z2;                           // PHI''(t)/PHI(t)
   double Zd;                           // derivative in Newton Raphson iteration
   double rdm1;                         // r * d - 1
   double tr;                           // 1/t
   double log2t;                        // log2(t)
   double method;                       // 0 for z2'(t) method, 1 for z3(t) method
   int i;                               // color
   int iter;                            // count iterations

   rdm1 = rd - 1.;
   if ( t_from == 0 && rdm1 <= 1. ) return 0.; //no inflection point
   rho[0] = r*omega;
   rho[1] = r;
   xx[0] = x;
   xx[1] = n - x;
   t = 0.5 * ( t_from + t_to );
   for ( i=0; i<COLORS; i++ ) {         // calculate zeta coefficients
      zeta[i][1][1] = rho[i];
      zeta[i][1][2] = rho[i] * ( rho[i] - 1. );
      zeta[i][2][2] = rho[i] * rho[i];
      zeta[i][1][3] = zeta[i][1][2] * ( rho[i] - 2. );
      zeta[i][2][3] = zeta[i][1][2] * rho[i] * 3.;
      zeta[i][3][3] = zeta[i][2][2] * rho[i] * 2.;
   }
   iter = 0;

   do {
      t1 = t;
      tr = 1. / t;
      log2t = log( t )*( 1./M_LN2 );
      phi[1] = phi[2] = phi[3] = 0.;
      for ( i=0; i<COLORS; i++ ) {       // calculate first 3 derivatives of phi(t)
         q1 = pow2_1( rho[i]*log2t,&q );
         q /= q1;
         phi[1] -= xx[i] * zeta[i][1][1] * q;
         phi[2] -= xx[i] * q * ( zeta[i][1][2] + q * zeta[i][2][2] );
         phi[3] -= xx[i] * q * ( zeta[i][1][3] + q * ( zeta[i][2][3] + q * zeta[i][3][3] ) );
      }
      phi[1] += rdm1;
      phi[2] -= rdm1;
      phi[3] += 2. * rdm1;
      phi[1] *= tr;
      phi[2] *= tr * tr;
      phi[3] *= tr * tr * tr;
      method = ( iter & 2 ) >> 1;        // alternate between the two methods
      Z2 = phi[1]*phi[1] + phi[2];
      Zd = method*phi[1]*phi[1]*phi[1] + ( 2.+method )*phi[1]*phi[2] + phi[3];

      if ( t < 0.5 ) {
         if ( Z2 > 0 ) {
            t_from = t;
         } else {
            t_to = t;
         }
         if ( Zd >= 0 ) {
            // use binary search if Newton-Raphson iteration makes problems
            t = ( t_from ? 0.5 : 0.2 ) * ( t_from + t_to );
         } else {
            // Newton-Raphson iteration
            t -= Z2 / Zd;
         }
      } else {
         if ( Z2 < 0 ) {
            t_from = t;
         } else {
            t_to = t;
         }
         if ( Zd <= 0 ) {
            // use binary search if Newton-Raphson iteration makes problems
            t = 0.5 * ( t_from + t_to );
         } else {
            // Newton-Raphson iteration
            t -= Z2 / Zd;
         }
      }
      if ( t >= t_to ) t = ( t1 + t_to ) * 0.5;
      if ( t <= t_from ) t = ( t1 + t_from ) * 0.5;
      if ( ++iter > 20 ) FatalError( "Search for inflection point failed in function search_inflect" );
   } while ( fabs( t - t1 ) > 1E-5 );
   return t;
}


double NonCentralHypergeometricC::probability( long int x_ ) {
   // calculate probability function. choosing best method
   x = x_;
   if ( x < xmin || x > xmax ) return 0.;
   if ( xmin == xmax ) return 1.;
   if ( omega == 1. ) { // hypergeometric
      return exp( lnbico() + LnFac( n ) + LnFac( N-n ) - LnFac( N ) );
   }

   long int x2 = n - x;
   long int x0 = x < x2 ? x : x2;
   int em = ( x == m || x2 == N-m );

   if ( x0 == 0 && n > 500 ) {
      return binoexpand();
   }

   if ( double( n )*x0 < 1000 || double( n )*x0 < 10000 && ( N > 1000.*n || em ) ) {
      return recursive();
   }

   if ( x0 <= 1 && N-n < 3 ) {
      return binoexpand();
   }

   findpars();
   if ( w < 0.04 && E < 10 && ( !em || w > 0.004 ) ) {
      return laplace();
   }

   return integrate();
}


int NonCentralHypergeometricC::MakeTable( double * table, long int MaxLength, long int * start, long int * end ) {
   // Makes a table of the noncentral hypergeometric probability distribution.
   // Returns true if table is long enough. Otherwise it fills the table
   // with as many correct values as possible and returns false.
   // The tails are cut off where the values are < 0.01*accuracy, so that
   // *start may be > xmin and *end may be < xmax.
   // The first and last x value represented in the table are returned in
   // *start and *end. The resulting probability values are returned in the
   // first (*end - *start + 1) positions of table. Any unused part of table
   // is overwritten with garbage. It is the responsibility of the calling
   // program to allocate sufficient space for the table.

   const int BUFSIZE = 2048;            // length of internal table
   double accuracya;                    // absolute accuracy
   double p[BUFSIZE+2];                 // probabilities
   double * p1, * p2;                   // offset into p
   double mxo;                          // (m-x)*omega
   double Nmnx;                         // N-m-nu+x
   double y, y1;                        // probability. Save old p[x] before it is overwritten
   double d1, d2, dcom;                 // divisors in probability formula
   long int xi, nu;                     // xi, nu = recursion values of x, n
   long int x1, x2;                     // lowest and highest x or xi
   long int i1, i2;                     // index into table
   double area;                         // estimate of area needed for recursion method

   area = double( n )*( xmax - xmin + 1 );
   if ( area < 4000 || area < 10000 && N > 1000.*n ) {
      // use recursion method
      accuracya = 0.005 * accuracy;      // absolute accuracy
      p1 = p2 = p + 1;                   // make space for p1[-1]
      p1[-1] = 0.;
      p1[0]  = 1.;         // initialize for recursion
      x1 = x2 = 0;
      for ( nu=1; nu<=n; nu++ ) {
         if ( n - nu < xmin - x1 || p1[x1] < accuracya ) {
            x1++;                          // increase lower limit when breakpoint passed or probability negligible
            p2--;
         }                         // compensate buffer offset in order to reduce storage space
         if ( x2 < xmax && p1[x2] >= accuracya ) {
            x2++;
            y1 = 0.;
         }               // increase upper limit until x has been reached
         else {
            y1 = p1[x2];
         }
         if ( p2 + x2 - p > BUFSIZE || x1 > x2 ) {
            goto ONE_BY_ONE;
         }              // Error: table length exceeded. Use other method

         mxo = ( m-x2 )*omega;
         Nmnx = N-m-nu+x2+1;
         for ( xi = x2; xi >= x1; xi-- ) { // backwards loop
            d2 = mxo + Nmnx;
            mxo += omega;
            Nmnx--;
            d1 = mxo + Nmnx;
            dcom = 1. / ( d1 * d2 );       // save a division by making common divisor
            y  = p1[xi-1]*mxo*d2*dcom + y1*( Nmnx+1 )*d1*dcom;
            y1 = p1[xi-1];                 // (warning: pointer alias, can't swap instruction order)
            p2[xi] = y;
         }
         p1 = p2;
      }

      // return results
      i1 = i2 = x2 - x1 + 1;             // desired table length
      if ( i2 > MaxLength ) i2 = MaxLength; // limit table length
      *start = x1;
      *end = x1 + i2 - 1;
      if ( i2>0 ) memcpy( table, p+1, i2*sizeof( table[0] ) ); // copy into external table
      return i1==i2;
   }                    // true if table size not reduced

   else {
      // area too big for recursion method
      // Calculate values one by one
ONE_BY_ONE:
      accuracya = 0.01 * accuracy;       // absolute accuracy

      // start to fill table from the end and down. start with x = floor(mean)
      x2 = ( long int )mean();
      x1 = x2 + 1;
      i1 = MaxLength;
      while ( x1 > xmin ) {              // loop for left tail
         x1--;
         i1--;
         y = probability( x1 );
         table[i1] = y;
         if ( y < accuracya ) break;
         if ( i1 == 0 ) break;
      }
      *start = x1;
      i2 = x2-x1+1;
      if ( i1>0 && i2>0 ) { // move numbers down to beginning of table
         memcpy( table, table+i1, i2*sizeof( table[0] ) );
      }
      // fill rest of table from mean and up
      i2--;
      while ( x2 < xmax ) {              // loop for right tail
         if ( i2 == MaxLength-1 ) {
            *end = x2;
            return 0;
         }          // table full
         x2++;
         i2++;
         y = probability( x2 );
         table[i2] = y;
         if ( y < accuracya ) break;
      }
      *end = x2;
      return 1;
   }
}


/***********************************************************************
         calculation methods in class MultiNonCentralHypergeometricC
***********************************************************************/

MultiNonCentralHypergeometricC::MultiNonCentralHypergeometricC( long int n_, long int * m_, double * odds_, int colors_, float accuracy_ ) {
   // constructor
   long int N1;
   int i;
   n = n_;
   m = m_;
   omega = odds_;
   colors = colors_;
   accuracy = accuracy_;
   r = 1.;
   for ( N=N1=0, i=0; i < colors; i++ ) {
      if ( m[i] < 0 || omega[i] < 0 ) FatalError( "Parameter negative in constructor for MultiNonCentralHypergeometricC" );
      N += m[i];
      if ( omega[i] ) N1 += m[i];
   }
   if ( N < n ) FatalError( "Not enough items in constructor for MultiNonCentralHypergeometricC" );
   if ( N1< n ) FatalError( "Not enough items with nonzero weight in constructor for MultiNonCentralHypergeometricC" );
}


void MultiNonCentralHypergeometricC::mean( double * mu ) {
   // calculate approximate mean of multivariate noncentral hypergeometric
   // distribution. Result is returned in mu[0..colors-1]
   double omeg[MAXCOLORS];              // scaled weights
   double omr;                          // reciprocal mean weight
   double t, t1;                        // independent variable in iteration
   double To, To1;                      // exp(t*omega[i]), 1-exp(t*omega[i])
   double H;                            // function to find root of
   double HD;                           // derivative of H
   int i;                               // color index
   int iter;                            // number of iterations

   // calculate mean weight
   for ( omr=0., i=0; i < colors; i++ ) omr += omega[i]*m[i];
   omr = N / omr;
   // scale weights to make mean = 1
   for ( i = 0; i < colors; i++ ) omeg[i] = omega[i] * omr;
   // Newton Raphson iteration
   iter = 0;
   t = -1.;                  // first guess
   do {
      t1 = t;
      H = HD = 0.;
      // calculate H and HD
      for ( i = 0; i < colors; i++ ) {
         if ( omeg[i] != 0. ) {
            To1 = pow2_1( t * ( 1./M_LN2 ) * omeg[i], &To );
            H += m[i] * To1;
            HD -= m[i] * omeg[i] * To;
         }
      }
      t -= ( H-n ) / HD;
      if ( t >= 0 ) t = 0.5 * t1;
      if ( ++iter > 20 ) {
         FatalError( "Search for mean failed in MultiNonCentralHypergeometricC::mean" );
      }
   } while ( fabs( H - n ) > 1E-3 );
   // finished iteration. Get all mu[i]
   for ( i = 0; i < colors; i++ ) {
      if ( omeg[i] != 0. ) {
         To1 = pow2_1( t * ( 1./M_LN2 ) * omeg[i] );
         mu[i] = m[i] * To1;
      } else {
         mu[i] = 0.;
      }
   }
}


// implementations of different calculation methods
double MultiNonCentralHypergeometricC::binoexpand( void ) {
   // binomial expansion of integrand
   // only implemented for x[i] = 0 for all but one i
   int i, j, k;
   double W = 0.;                       // total weight
   for ( i=j=k=0; i<colors; i++ ) {
      W += omega[i] * m[i];
      if ( x[i] ) {
         j=i;
         k++;
      }
   }                      // find the nonzero x[i]
   if ( k > 1 ) FatalError( "More than one x[i] nonzero in MultiNonCentralHypergeometricC::binoexpand" );
   return exp( FallingFactorial( m[j],n ) - FallingFactorial( W/omega[j],n ) );
}


double MultiNonCentralHypergeometricC::lnbico( void ) {
   // natural log of binomial coefficients
   bico = 0.;
   int i;
   for ( i=0; i<colors; i++ ) {
      if ( x[i] < m[i] && omega[i] ) {
         bico += LnFac( m[i] ) - LnFac( x[i] ) - LnFac( m[i]-x[i] );
      }
   }
   return bico;
}


void MultiNonCentralHypergeometricC::findpars( void ) {
   // calculate r, w, E
   // calculate d, E, r, w

   // find r to center peak of integrand at 0.5
   double dd;                           // scaled d
   double dr;                           // 1/d

   double z, zd, rr, lastr, rrc, rt, r2, r21, a, b, ro, k1;
   double omax;                         // highest omega
   double omaxr;                        // 1/omax
   double omeg[MAXCOLORS];              // scaled weights
   int i, j = 0;

   // find highest omega
   for ( omax=0., i=0; i < colors; i++ ) {
      if ( omega[i] > omax ) omax = omega[i];
   }
   omaxr = 1. / omax;
   dd = E = 0.;
   for ( i = 0; i < colors; i++ ) {
      // scale weights to make max = 1
      omeg[i] = omega[i] * omaxr;
      // calculate d and E
      dd += omeg[i] * ( m[i]-x[i] );
      E  += omeg[i] * m[i];
   }
   dr = 1. / dd;
   E *= dr;
   rr = r * omax;
   if ( rr <= dr ) rr = 1.2 * dr; // initial guess
   // Newton-Raphson iteration to find r
   do {
      lastr = rr;
      rrc = 1. / rr;
      z = dd - rrc;                      // z(r)
      zd = rrc * rrc;                    // z'(r)
      for ( i=0; i<colors; i++ ) {
         rt = rr * omeg[i];
         if ( rt < 100. && rt > 0. ) {    // avoid overflow and division by 0
            r21 = pow2_1( rt, &r2 );       // r2=2^r, r21=1.-2^r
            a = omeg[i] / r21;             // omegai/(1.-2^r)
            b = x[i] * a;                  // x*omegai/(1.-2^r)
            z  += b;
            zd += b * a * r2 * M_LN2;
         }
      }
      if ( zd == 0 ) FatalError( "can't find r" );
      rr -= z / zd;                      // next r
      if ( rr <= dr ) rr = lastr * 0.125 + dr * 0.875;
      if ( ++j == 70 ) FatalError( "convergence problem searching for r" );
   } while ( fabs( rr-lastr ) > rr * 1.E-5 );
   rd = rr * dd;
   r = rr * omaxr;

   // find peak width
   phi2d = 0.;
   for ( i=0; i<colors; i++ ) {
      ro = rr * omeg[i];
      if ( ro < 300 && ro > 0. ) {       // avoid overflow and division by 0
         k1 = pow2_1( ro );
         k1 = -1. / k1;
         k1 = omeg[i] * omeg[i] * ( k1 + k1*k1 );
      } else k1 = 0.;
      phi2d += x[i] * k1;
   }
   phi2d *= -4. * rr * rr;
   if ( phi2d > 0. ) FatalError( "peak width undefined in function findpars" );
   w = 1./sqrt( -phi2d );
}


double MultiNonCentralHypergeometricC::laplace( void ) {
   // Laplace's method with narrow integration interval,
   // using error function residue table, defined in erfres.cpp
   // NOTE: This function can only be used when the integrand peak is narrow.
   // findpars() must be called before this function.

   const int MAXDEG = 40;        // arraysize
   int degree;                   // max expansion degree
   double accur;                 // stop expansion when terms below this threshold
   double f0;                    // factor outside integral
   double rho[MAXCOLORS];        // r*omegai
   double qi;                    // 2^(-rho)
   double qi1;                   // 1-qi
   double qq[MAXCOLORS];         // qi / qi1
   double eta[MAXCOLORS+1][MAXDEG+1]; // eta coefficients
   double phideri[MAXDEG+1];     // derivatives of phi
   double PSIderi[MAXDEG+1];     // derivatives of PSI
   double epsilon;               // factor for integration interval
   double epsilonmax;            // max possible epsilon
   double * erfres;              // pointer to table of error function expansion residues
   int erfreslen;                // length of table

   // variables in asymptotic summation
   static const double sqrtpi = 1.772453850905515993; // sqrt(pi)
   static const double sqrt8  = 2.828427124746190098; // sqrt(8)
   double qqpow;                 // qq^j
   double pow2k;                 // 2^k
   double mfac;                  // (-1)^(k-1)*(k-1)!
   double kfac;                  // k!
   double bino;                  // binomial coefficient
   double wr;                    // 1/w, w = integrand peak width
   double vr;                    // 1/v, v = integration interval
   double v2m2;                  // (2*v)^(-2)
   double v2mk1;                 // (2*v)^(-k-1)
   double gammal2;               // Gamma((k+1)/2);
   double s;                     // summation term
   double sum;                   // Taylor sum

   int i;                        // loop counter for color
   int j;                        // loop counter for derivative
   int k;                        // loop counter for expansion degree
   int ll;                       // k/2
   int converg=0;                // number of consequtive terms below accuracy

   // initialize
   for ( k=0; k<=2; k++ )  phideri[k] = PSIderi[k] = 0;

   // find rho[i], qq[i], first eta coefficients, and zero'th derivative of phi
   for ( i=0; i<colors; i++ ) {
      rho[i] = r * omega[i];
      if ( rho[i] == 0. ) continue;
      if ( rho[i] > 40. ) {
         qi=0.;
         qi1 = 1.;
      }  // avoid underflow
      else {
         qi1 = pow2_1( -rho[i], &qi );
      }   // qi=2^(-rho), qi1=1.-2^(-rho)
      qq[i] = qi / qi1;    // 2^(-r*omegai)/(1.-2^(-r*omegai))
      // peak = zero'th derivative
      phideri[0] += x[i] * log1mx( qi, qi1 );
      // eta coefficients
      eta[i][0] = 0.;
      eta[i][1] = eta[i][2] = rho[i]*rho[i];
   }

   // d, r, and w must be calculated by findpars()
   // zero'th derivative
   phideri[0] -= ( rd - 1. ) * M_LN2;
   // scaled factor outside integral
   f0 = rd * exp( phideri[0] + lnbico() );
   // calculate narrowed integration interval
   vr = sqrt8 * w;
   wr = 1. / w;
   phideri[2] = phi2d;
   epsilonmax = 0.5 * wr;
   if ( accuracy < 1.E-10 && epsilonmax >= 8. ) {
      epsilon = 8.;
      erfres = erfres80;
      erfreslen = sizeof( erfres80 )/sizeof( erfres80[0] );
   } else if ( accuracy < 1.E-8 && epsilonmax >= 7. ) {
      epsilon = 7.;
      erfres = erfres70;
      erfreslen = sizeof( erfres70 )/sizeof( erfres70[0] );
   } else if ( accuracy < 1.E-6 && epsilonmax >= 6. ) {
      epsilon = 6.;
      erfres = erfres60;
      erfreslen = sizeof( erfres60 )/sizeof( erfres60[0] );
   } else {
      epsilon = 5.;
      erfres = erfres50;
      erfreslen = sizeof( erfres50 )/sizeof( erfres50[0] );
   }
   if ( epsilon > epsilonmax ) FatalError( "Laplace failed. peak width too high" );
   degree = MAXDEG;
   if ( degree >= erfreslen*2 ) degree = erfreslen*2-2;

   // variables used in summation
   v2m2 = 0.25*vr*vr;        // (2*v)^(-2)

   // set up for starting loop at k=3
   PSIderi[0] = 1.;
   pow2k = 8.;
   mfac = 2.;
   kfac = 2.;
   sum = 0.5 * vr * sqrtpi * erfres[0];
   v2mk1 = 0.5*vr*v2m2*v2m2;
   gammal2 = sqrtpi * 0.75;
   accur = accuracy * sum;

   // summation loop
   for ( k=3; k <= degree; k++ ) {
      phideri[k] = 0.;

      // loop for all colors
      for ( i=0; i<colors; i++ ) {
         if ( rho[i] == 0. ) continue;
         eta[i][k] = 0.;
         // backward loop for all powers
         for ( j=k; j>0; j-- ) {
            // find coefficients recursively from previous coefficients
            eta[i][j]  =  eta[i][j]*( j*rho[i]-( k-2 ) ) +  eta[i][j-1]*rho[i]*( j-1 );
         }
         qqpow = 1.;
         // forward loop for all powers
         for ( j=1; j<=k; j++ ) {
            qqpow *= qq[i];   // qq^j
            // contribution to derivative
            phideri[k] += x[i] * eta[i][j] * qqpow;
         }
      }

      // finish calculation of derivatives
      phideri[k] = -pow2k * phideri[k] + 2*( 1-k )*phideri[k-1];

      pow2k *= 2.;    // 2^k
      mfac *= -k;     // (-1)^(k-1)*(k-1)!

      // loop to calculate derivatives of PSI from derivatives of psi.
      // terms # 0, 1, 2, k-2, and k-1 are zero and not included in loop.
      // The j'th derivatives of psi are identical to the derivatives of phi for j>2, and
      // zero for j=1,2. Hence we are using phideri[j] for j>2 here.
      PSIderi[k] = phideri[k];  // this is term # k
      bino = 0.5*( k-1 )*( k-2 ); // binomial coefficient for term # 3
      for ( j=3; j < k-2; j++ ) { // loop for remaining nonzero terms (if k>5)
         PSIderi[k] += PSIderi[k-j] * phideri[j] * bino;
         bino *= double( k-j )/double( j );
      }

      kfac *= k;

      if ( ( k & 1 ) == 0 ) { // only for even k
         ll = k/2;
         s = ( PSIderi[k]/kfac )*v2mk1* gammal2 * erfres[ll];
         sum += s;

         // check for convergence of Taylor expansion
         if ( fabs( s ) < accur ) converg++;
         else converg = 0;
         if ( converg > 1 ) break;

         // update recursive expressions
         v2mk1 *= v2m2;
         gammal2 *= k/2 + 0.5;
      }
   }

   // multiply by terms outside integral
   return f0 * sum;
}


double MultiNonCentralHypergeometricC::integrate( void ) {
   // non-central hypergeometric distribution function
   // calculation by numerical integration with variable-length steps
   // NOTE: findpars() must be called before this function.
   double s;                            // result of integration step
   double sum;                          // integral
   double ta, tb;                       // subinterval for integration step

   lnbico();                            // compute log of binomial coefficients

   // choose method:
   if ( w < 0.02 ) {
      // normal method. Step length determined by peak width w
      double delta, s1;
      s1 = accuracy < 1E-9 ? 0.5 : 1.;
      delta = s1 * w;                    // integration steplength
      ta = 0.5 + 0.5 * delta;
      sum = integrate_step( 1.-ta, ta ); // first integration step around center peak
      do {
         tb = ta + delta;
         if ( tb > 1. ) tb = 1.;
         s  = integrate_step( ta, tb );   // integration step to the right of peak
         s += integrate_step( 1.-tb,1.-ta ); // integration step to the left of peak
         sum += s;
         if ( s < accuracy * sum ) break; // stop before interval finished if accuracy reached
         ta = tb;
         if ( tb > 0.5 + w ) delta *= 2.;
      }  // increase step length far from peak
      while ( tb < 1. );
   }

   else {
      // difficult situation. Step length determined by inflection points
      double t1, t2, tinf, delta, delta1;
      sum = 0.;
      // do left and right half of integration interval separately:
      for ( t1=0., t2=0.5; t1 < 1.; t1+=0.5, t2+=0.5 ) {
         // integrate from 0 to 0.5 or from 0.5 to 1
         tinf = search_inflect( t1, t2 ); // find inflection point
         delta = tinf - t1;
         if ( delta > t2 - tinf ) delta = t2 - tinf; // distance to nearest endpoint
         delta *= 1./7.;                  // 1/7 will give 3 steps to nearest endpoint
         if ( delta < 1E-4 ) delta = 1E-4;
         delta1 = delta;
         // integrate from tinf forwards to t2
         ta = tinf;
         do {
            tb = ta + delta1;
            if ( tb > t2 - 0.25*delta1 ) tb = t2; // last step of this subinterval
            s = integrate_step( ta, tb );  // integration step
            sum += s;
            delta1 *= 2;                   // double steplength
            if ( s < sum * 1E-4 ) delta1 *= 8.; // large step when s small
            ta = tb;
         } while ( tb < t2 );
         if ( tinf ) {
            // integrate from tinf backwards to t1
            tb = tinf;
            do {
               ta = tb - delta;
               if ( ta < t1 + 0.25*delta ) ta = t1; // last step of this subinterval
               s = integrate_step( ta, tb ); // integration step
               sum += s;
               delta *= 2;                  // double steplength
               if ( s < sum * 1E-4 ) delta *= 8.; // large step when s small
               tb = ta;
            } while ( ta > t1 );
         }
      }
   }

   return sum * rd;
}


double MultiNonCentralHypergeometricC::integrate_step( double ta, double tb ) {
   // integration subprocedure used by integrate()
   // makes one integration step from ta to tb using Gauss-Legendre method.
   // result is scaled by multiplication with exp(bico)
   double ab, delta, tau, ltau, y, sum, taur, rdm1;
   int i, j;

   // define constants for Gauss-Legendre integration with IPOINTS points
#define IPOINTS  8  // number of points in each integration step

#if IPOINTS == 3
   static const double xval[3]    = {-.774596669241,0,0.774596668241};
   static const double weights[3] = {.5555555555555555,.88888888888888888,.55555555555555};
#elif IPOINTS == 4
   static const double xval[4]    = {-0.861136311594,-0.339981043585,0.339981043585,0.861136311594},
                                    static const double weights[4] = {0.347854845137,0.652145154863,0.652145154863,0.347854845137};
#elif IPOINTS == 5
   static const double xval[5]    = {-0.906179845939,-0.538469310106,0,0.538469310106,0.906179845939};
   static const double weights[5] = {0.236926885056,0.478628670499,0.568888888889,0.478628670499,0.236926885056};
#elif IPOINTS == 6
   static const double xval[6]    = {-0.932469514203,-0.661209386466,-0.238619186083,0.238619186083,0.661209386466,0.932469514203};
   static const double weights[6] = {0.171324492379,0.360761573048,0.467913934573,0.467913934573,0.360761573048,0.171324492379};
#elif IPOINTS == 8
   static const double xval[8]    = {-0.960289856498,-0.796666477414,-0.525532409916,-0.183434642496,0.183434642496,0.525532409916,0.796666477414,0.960289856498};
   static const double weights[8] = {0.10122853629,0.222381034453,0.313706645878,0.362683783378,0.362683783378,0.313706645878,0.222381034453,0.10122853629};
#elif IPOINTS == 12
   static const double xval[12]   = {-0.981560634247,-0.90411725637,-0.769902674194,-0.587317954287,-0.367831498998,-0.125233408511,0.125233408511,0.367831498998,0.587317954287,0.769902674194,0.90411725637,0.981560634247};
   static const double weights[12]= {0.0471753363866,0.106939325995,0.160078328543,0.203167426723,0.233492536538,0.249147045813,0.249147045813,0.233492536538,0.203167426723,0.160078328543,0.106939325995,0.0471753363866};
#elif IPOINTS == 16
   static const double xval[16]   = {-0.989400934992,-0.944575023073,-0.865631202388,-0.755404408355,-0.617876244403,-0.458016777657,-0.281603550779,-0.0950125098376,0.0950125098376,0.281603550779,0.458016777657,0.617876244403,0.755404408355,0.865631202388,0.944575023073,0.989400934992};
   static const double weights[16]= {0.027152459411,0.0622535239372,0.0951585116838,0.124628971256,0.149595988817,0.169156519395,0.182603415045,0.189450610455,0.189450610455,0.182603415045,0.169156519395,0.149595988817,0.124628971256,0.0951585116838,0.0622535239372,0.027152459411};
#else
#error // IPOINTS must be a value for which the tables are defined
#endif

   delta = 0.5 * ( tb - ta );
   ab = 0.5 * ( ta + tb );
   rdm1 = rd - 1.;
   sum = 0;

   for ( j = 0; j < IPOINTS; j++ ) {
      tau = ab + delta * xval[j];
      ltau = log( tau );
      taur = r * ltau;
      y = 0.;
      for ( i = 0; i < colors; i++ ) {
         // possible loss of precision due to subtraction here:
         if ( omega[i] ) {
            y += pow1pow( taur*omega[i],x[i] );
         }
      } // ln((1-e^taur*omegai)^xi)
      y += rdm1*ltau + bico;
      if ( y > -50. ) sum += weights[j] * exp( y );
   }
   return delta * sum;
}


double MultiNonCentralHypergeometricC::search_inflect( double t_from, double t_to ) {
   // search for an inflection point of the integrand PHI(t) in the interval
   // t_from < t < t_to
   double t, t1;                        // independent variable
   double rho[MAXCOLORS];               // r*omega[i]
   double q;                            // t^rho[i] / (1-t^rho[i])
   double q1;                           // 1-t^rho[i]
   double zeta[MAXCOLORS][4][4];        // zeta[i,j,k] coefficients
   double phi[4];                       // derivatives of phi(t) = log PHI(t)
   double Z2;                           // PHI''(t)/PHI(t)
   double Zd;                           // derivative in Newton Raphson iteration
   double rdm1;                         // r * d - 1
   double tr;                           // 1/t
   double log2t;                        // log2(t)
   double method;                       // 0 for z2'(t) method, 1 for z3(t) method
   int i;                               // color
   int iter;                            // count iterations

   rdm1 = rd - 1.;
   if ( t_from == 0 && rdm1 <= 1. ) return 0.; //no inflection point
   t = 0.5 * ( t_from + t_to );
   for ( i=0; i<colors; i++ ) {         // calculate zeta coefficients
      rho[i] = r * omega[i];
      zeta[i][1][1] = rho[i];
      zeta[i][1][2] = rho[i] * ( rho[i] - 1. );
      zeta[i][2][2] = rho[i] * rho[i];
      zeta[i][1][3] = zeta[i][1][2] * ( rho[i] - 2. );
      zeta[i][2][3] = zeta[i][1][2] * rho[i] * 3.;
      zeta[i][3][3] = zeta[i][2][2] * rho[i] * 2.;
   }
   iter = 0;

   do {
      t1 = t;
      tr = 1. / t;
      log2t = log( t )*( 1./M_LN2 );
      phi[1] = phi[2] = phi[3] = 0.;
      for ( i=0; i<colors; i++ ) {       // calculate first 3 derivatives of phi(t)
         if ( rho[i] == 0. ) continue;
         q1 = pow2_1( rho[i]*log2t,&q );
         q /= q1;
         phi[1] -= x[i] * zeta[i][1][1] * q;
         phi[2] -= x[i] * q * ( zeta[i][1][2] + q * zeta[i][2][2] );
         phi[3] -= x[i] * q * ( zeta[i][1][3] + q * ( zeta[i][2][3] + q * zeta[i][3][3] ) );
      }
      phi[1] += rdm1;
      phi[2] -= rdm1;
      phi[3] += 2. * rdm1;
      phi[1] *= tr;
      phi[2] *= tr * tr;
      phi[3] *= tr * tr * tr;
      method = ( iter & 2 ) >> 1;        // alternate between the two methods
      Z2 = phi[1]*phi[1] + phi[2];
      Zd = method*phi[1]*phi[1]*phi[1] + ( 2.+method )*phi[1]*phi[2] + phi[3];

      if ( t < 0.5 ) {
         if ( Z2 > 0 ) {
            t_from = t;
         } else {
            t_to = t;
         }
         if ( Zd >= 0 ) {
            // use binary search if Newton-Raphson iteration makes problems
            t = ( t_from ? 0.5 : 0.2 ) * ( t_from + t_to );
         } else {
            // Newton-Raphson iteration
            t -= Z2 / Zd;
         }
      } else {
         if ( Z2 < 0 ) {
            t_from = t;
         } else {
            t_to = t;
         }
         if ( Zd <= 0 ) {
            // use binary search if Newton-Raphson iteration makes problems
            t = 0.5 * ( t_from + t_to );
         } else {
            // Newton-Raphson iteration
            t -= Z2 / Zd;
         }
      }
      if ( t >= t_to ) t = ( t1 + t_to ) * 0.5;
      if ( t <= t_from ) t = ( t1 + t_from ) * 0.5;
      if ( ++iter > 20 ) FatalError( "Search for inflection point failed in function search_inflect" );
   } while ( fabs( t - t1 ) > 1E-5 );
   return t;
}


double MultiNonCentralHypergeometricC::probability( long int * x_ ) {
   // calculate probability function. choosing best method
   int i, j, em;
   long int xsum;
   x = x_;

   for ( xsum=i=0; i<colors; i++ )  xsum += x[i];
   if ( xsum != n )
      FatalError( "sum of x values not equal to n in function MultiNonCentralHypergeometricC::probability" );

   if ( colors < 3 ) {
      if ( colors <= 0 ) return 1.;
      if ( colors == 1 ) return x[0] == m[0];
      // colors = 2
      if ( omega[1] == 0. ) return x[0] == m[0];
      return NonCentralHypergeometricC( n,m[0],N,omega[0]/omega[1],accuracy ).probability( x[0] );
   }

   for ( i=j=em=0; i<colors; i++ ) {
      if ( x[i] > m[i] ) return 0.;
      if ( omega[i]==0. && x[i] ) return 0.;
      if ( x[i] > 0 ) j++;
      if ( x[i] == m[i] ) em++;
   }

   if ( n == 0 || n == N ) return 1.;

   if ( j == 1 )
      return binoexpand();

   findpars();
   if ( w < 0.04 && E < 10 && ( !em || w > 0.004 ) ) {
      return laplace();
   }

   return integrate();
}


/***********************************************************************
         Methods for MultiNonCentralHypergeometricMoments
***********************************************************************/

double MultiNonCentralHypergeometricMoments::moments( double * mu, double * variance, long int * combinations ) {
   // calculates mean and variance of the noncentral hypergeometric
   // distribution by calculating all combinations of x-values.
   // Return value = sum of all probabilities. The deviation of this value
   // from 1 is a measure of the accuracy.
   // Returns the mean to mean[0...colors-1]
   // Returns the variance to variance[0...colors-1]
   double sumf;               // sum of all f(x) values
   long int msum;             // temporary sum
   int i;                     // loop counter

   // get approximate mean
   mean( sx );
   // round mean to integers
   for ( i=0; i < colors; i++ ) {
      xm[i] = ( long int )( sx[i]+0.4999999 );
   }

   // set up for recursive loops
   for ( i=colors-1, msum=0; i >= 0; i-- ) {
      remaining[i] = msum;
      msum += m[i];
   }
   for ( i=0; i<colors; i++ )  sx[i] = sxx[i] = 0.;
   sn = 0;

   // recursive loops to calculate sums
   sumf = loop( n, 0 );

   // calculate mean and variance
   for ( i = 0; i < colors; i++ ) {
      mu[i] = sx[i]/sumf;
      variance[i] = sxx[i]/sumf - sx[i]*sx[i]/( sumf*sumf );
   }

   // return combinations and sum
   if ( combinations ) *combinations = sn;
   return sumf;
}


double MultiNonCentralHypergeometricMoments::loop( long int n, int c ) {
   // recursive function to loop through all combinations of x-values.
   // used by moments()
   long int x, x0;                      // x of color c
   long int xmin, xmax;                 // min and max of x[c]
   double s1, s2, sum = 0.;             // sum of f(x) values
   int i;                               // loop counter

   if ( c < colors-1 ) {
      // not the last color
      // calculate min and max of x[c] for given x[0]..x[c-1]
      xmin = n - remaining[c];
      if ( xmin < 0 ) xmin = 0;
      xmax = m[c];
      if ( xmax > n ) xmax = n;
      x0 = xm[c];
      if ( x0 < xmin ) x0 = xmin;
      if ( x0 > xmax ) x0 = xmax;
      // loop for all x[c] from mean and up
      for ( x = x0, s2 = 0.; x <= xmax; x++ ) {
         xi[c] = x;
         sum += s1 = loop( n-x, c+1 ); // recursive loop for remaining colors
         if ( s1 < accuracy && s1 < s2 ) break; // stop when values become negligible
         s2 = s1;
      }
      // loop for all x[c] from mean and down
      for ( x = x0-1; x >= xmin; x-- ) {
         xi[c] = x;
         sum += s1 = loop( n-x, c+1 ); // recursive loop for remaining colors
         if ( s1 < accuracy && s1 < s2 ) break; // stop when values become negligible
         s2 = s1;
      }
   } else {
      // last color
      xi[c] = n;
      s1 = probability( xi );
      for ( i=0; i < colors; i++ ) {
         sx[i]  += s1 * xi[i];
         sxx[i] += s1 * xi[i] * xi[i];
      }
      sn++;
      sum = s1;
   }
   return sum;
}
